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The aeroelastic and aerothermoelastic behavior of three-dimensional configurations in 
hypersonic flow regime are studied. The aeroelastic behavior of a low aspect ratio wing, 
representative of a fin or control surface on a generic hypersonic vehicle, is examined using 
third order piston theory, Euler and Navier-Stokes aerodynamics. The sensitivity of the 
aeroelastic behavior generated using Euler and Navier-Stokes aerodynamics to parameters 
governing temporal accuracy is also examined. Also, a refined aerothermoelastic model, 
which incorporates the heat transfer between the fluid and structure using CFD generated 
aerodynamic heating, is used to examine the aerothermoelastic behavior of the low as- 
pect ratio wing in the hypersonic regime. Finally, the hypersonic aeroelastic behavior of a 
generic hypersonic vehicle with a lifting-body type fuselage and canted fins is studied using 
piston theory and Euler aerodynamics for the range of 2.5 < M < 28, at altitudes ranging 
from 10,000 feet to 80,000 feet. This analysis includes a study on optimal mesh selection for 
use with Euler aerodynamics. In addition to the aeroelastic and aerothermoelastic results 
presented, three time domain flutter identification techniques are compared, namely the 
moving block approach, the least squares curve fitting method, and a system identification 
technique using an Auto-Regressive model of the aeroelastic system. In general, the three 
methods agree well. The system identification technique, however, provided quick damp- 
ing and frequency estimations with minimal response record length, and therefore offers 
significant reductions in computational cost. In the present case, the computational cost 
was reduced by 75%. The aeroelastic and aerothermoelastic results presented illustrate 
the applicability of the CFL3D code for the hypersonic flight regime. 
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Nondimensional offset between the elastic axis and the midchord, positive for elastic axis 
behind midchord 
Speed of sound 

Constant in sinusoidal representation of transient modal amplitude 
Estimated aeroelastic system matrix 

Coefficients used for damping and frequency identification 
Semi-chord 

Reference length, chord length of double- wedge airfoil 
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CoefRcients of lift, moment about the elastic axis, and drag 

Coefficient of pressure 

Estimated aeroelastic system matrix 

Specific heat of the wall 

CFL3D input parameter regulating “pseudo time step size 

Squared error between curve fit and actual data 

Exponential function for curve fit 

Airfoil vertical displacement at Elastic Axis 

States in state-space representation of AR model 

Heat transfer coefficient 

Mass moment of inertia about the Elastic Axis 
Discrete time 

Spring constants in pitch and plunge respectively; K a = I a u}^,Kh = mwj 

Reduced frequency 

Free stream Mach number 

Mass per unit span 

Flutter Mach number 

Generalized mass and stiffness matrices of the structure 
Number of points in response signal 
Normal vector 

Number of points in moving blocknumber 

Number of modes used 

Pressure 

Free-stream pressure 
Dynamic pressure 

Generalized force vector for the structure 
Generalized force corresponding to mode i 
Modal amplitude of mode i 

Heat transfer rate due to aerodynamic heating, radiation, conduction, 

and stored enegry respectively 

Finite Fourier transform of transient of mode i 

Amplitude of Fourier spectrum at one time-step (moving-block function) 

Real part of eigenvalue 

Nondimensional radius of gyration 

Reynolds number 

Surface area of the structure 

Imaginary part of eigenvalue 

Static mass moment of wing section about elastic axis 

Temperature 

Time 

Time step 

Kinetic energy of the structure 

Sample time 

Airfoil half thickness 

kth point in discrete time 

Radiation equilibrium wall temperature 

Response signal length 

Wall temperature 

Adiabatic wall temperature 

Reference temperature 

Potential energy of the structure 

Free stream velocity 

Normal velocity of airfoil surfaces 

Displacement of the surface of the structure 

Spatial Coordinates 
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Nondimensional offset between the elastic axis and the cross-sectional center of gravity, 
positive for center of gravity behind elastic axis 
State matrix 

Position of structural surface 

Airfoil pitch displacement about the Elastic Axis 

Static angle of attack 

Skin thickness 

Nondimensional grid-spacing 

Input for ARMA model of aeroelastic system 

Emissivity 

Ratio of specific heats 
Thermal conductivity 
Estimated matrix eigenvalue 
Mass ratio 

Natural frequencies of uncoupled pitch and plunge motions 

Frequency of mode i 

Natural frequency of mode i 

Modal matrix 

Mode shape for mode i 

Phase angle in radians of modal amplitude history in sinusoidal representation, for mode i 
Air density 
Wall density 

Stefan-Boltzmann constant 
Local inclination to the flow 
Time 

Damping ratio 

First and second derivatives with respect to time 


I. Introduction 

H ypersonic aeroelasticity and aerothermoelasticity was an active area of research in the late 1950’s and 
during the 1960’s as evident from Refs. 1-4. This research was instrumental in providing the basis for 
the aerothermoelastic design of the space shuttle. For a considerable time period there was only limited 
interest in this area until the advent of the National AeroSpace Plane (NASP). 

In recent years, renewed activity in hypersonic flight research has been stimulated by the need for a low 
cost, reusable launch vehicle (RLV) and the recently completed Hyper-X program at NASA. Other activities 
are motivated by the design requirements for unmanned hypersonic vehicles that meet the needs of the US 
Air Force. 

Vehicles in this category are based on a lifting body design. However, stringent minimum- weight re- 
quirements imply a degree of fuselage flexibility. Aerodynamic surfaces, needed for control, are also flexible. 
Furthermore, to meet the requirement of a flight profile that spans the Mach number range from 0 to 15, the 
vehicle must withstand severe aerodynamic heating. These factors combine to produce unusual aeroelastic 
problems that have received only limited attention in the past. Furthermore, it is important to emphasize 
that testing of aeroelastically scaled wind tunnel models, a conventional practice in subsonic and supersonic 
flow, is not feasible in the hypersonic regime. Thus, the role of aeroelastic simulations is critical for this 
flight regime. 

Previous studies in this area can be separated into several groups. The first group consists of studies 
focusing on panel flutter, which is a localized aeroelastic problem representing a small portion of the skin on 
the surface of the hypersonic vehicle. 5-10 The second group of studies in this area was motivated by a previous 
hypersonic vehicle, namely the NASP. 11-17 However, some of these studies dealt with the transonic regime, 
because it was perceived to be quite important. Spain et al. 12 carried out a flutter analysis of all-movable 
N ASP-like wings with slab and double- wedge airfoils, and found that using effective shapes for the airfoils 
obtained by adding the boundary layer displacement thickness to the airfoil thickness improved the overall 
agreement with experiments. Aerothermoelastic analyses of NASP-like vehicles found that aerodynamic 
heating altered the aeroelastic stability of the vehicle through the degradation of material properties and 
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introduction of thermal stresses. 15-17 

The third group of studies is restricted to recent papers that deal with the newer hypersonic configurations 
such as the X-33 or the X-34. Reference 18 considered the X-34 launch vehicle in free flight at M = 8.0. In 
Ref. 19, CFD-based flutter analysis was used for the aeroelastic analysis of the X-43 configuration, using 
system identification based order reduction of the aerodynamic degrees of freedom. Both the structure and 
the fluid were discretized using the finite element approach. It was shown that piston theory and AR.MA 
Euler calculations predicted somewhat similar results. In Ref. 20 the aeroelastic stability of a generic 
hypersonic vehicle resembling the X-33 was considered using piston theory. In a sequel to this study 21 
the piston theory results for the aeroelastic behavior were also compared to aeroelastic stability boundaries 
obtained from computational aeroelasticity. 

From the studies on previous hypersonic vehicles, 11,19 ’ 22 ’ 23 one can identify operating envelopes for each 
vehicle. A graphical representation of these operating conditions is shown in Fig. 1. This figure is a 
convenient illustration of typical operating conditions of generic hypersonic vehicles. 

The authors of this paper described an aeroelastic analysis capability for generic hypersonic vehicles in 
the Mach number range 0.5 < M < 15.0, using computational aeroelasticity. 24-26 The computational tool 
consisted of a combination of the CFL3D code 27 and a finite element model of both a generic hypersonic 
vehicle, and three-dimensional wing, utilizing MSC.NASTRAN. This paper continues to explore fundamental 
aspects of hypersonic aeroelasticity and aerothermoelasticity using computational tools and it focuses on 
three-dimensional low aspect ratio wing and complete generic hypersonic vehicle configurations. The specific 
objectives of the paper are: 

1. Compare the aeroelastic behavior predicted on a low aspect ratio wing using CFL3D, with both Euler 
and Navier-Stokes solutions for the unsteady aerodynamic loads. Examine the sensitivity of these 
computations to parameters relating to temporal accuracy. 

2. Study different time domain damping and frequency identification techniques in order to asses the 
properties and effectiveness of these methods. 

3. Obtain the stability boundaries for the low-aspect ratio wing using unsteady Euler and Navier-Stokes 
aerodynamics, and compare them with similar results based on third order piston theory. 

4. Develope a refined aerothermoelastic model that incorporates the effect of heat transfer between the 
structure and the fluid due to intense aerodynamic heating. 

5. Predict the aeroelastic behavior of a generic hypersonic vehicle using nonlinear piston theory and Euler 
aerodynamics over a broad range of Mach numbers and altitudes. Examine the sensitivity of the flutter 
prediction to grid resolution. 

These objectives will significantly enhance our understanding of hypersonic aeroelasticity as applied to 
reusable launch vehicles. 


II. Method of Solution 

The computational aeroelastic solutions in the present study are obtained using the CFL3D code. 2 ' The 
CFL3D code is used to generate both steady and unsteady air loads, and it also produces the aeroelastic 
transients and response solutions. The fluid/structure coupling in the code is accomplished using the free 
vibration modes of the vehicle. 

A. Euler/Navier-Stokes Solver in CFL3D 

The aeroelastic analysis in this study is carried out using the CFL3D code. The code uses an implicit, 
finite-volume algorithm based on upwind-biased spatial differencing to solve the time-dependent Euler and 
Reynolds-averaged Navier-Stokes equations. Multigrid and mesh-sequencing are available for convergence 
acceleration. The algorithm, which is based on a cell-centered scheme, uses upwind-differencing based on 
either flux-vector splitting or flux-difference splitting, and can sharply capture shock waves. For applications 
utilizing the thin-layer Navier-Stokes equations, different turbulence models are available. For time-accurate 
problems using a deforming mesh, an additional term accounting for the change in cell- volume is included 
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in the time-discretization of the governing equations. Since CFL3D is an implicit code using approximate 
factorization, linearization and factorization errors are introduced at every time-step. Hence, intermediate 
calculations referred to as “subiterations” are used to reduce these errors. Increasing these subiterations 
improves the accuracy of the simulation, albeit at increased computational cost. 


B. Aeroelastic Option in CFL3D 

The aeroelastic approach underlying the CFL3D code is similar to that described in Refs. 28 and 29. The 
equations are derived by assuming that the general motion w(x,y, z,t) of the structure is described by a 
finite modal series given by Eq. (1) below. The functions <pi(x,y,z) represent the free vibration modes of 
the vehicle which are calculated using a finite element approach. 

rim 

w(x,y,z,t) = '^2q i (t)(/>i(x,y,z) (1) 

*= l 


The aeroelastic equations of motion are obtained from Lagrange’s equations, 


_d ( OTe] _ 9Te + 9Ue 

dt V dqi ) dqi dqi 


i — 1,2,..., n m 


( 2 ) 


which 

where 


yield 

Mq + Kq = Q(q, q, q), q T = [<?i <?2 9nJ 
the elements of the generalized force vector are given by, 
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( 3 ) 

( 4 ) 


The aeroelastic equations are written in terms of a linear state-space equation (using a state vector of the 
form [... qi-\ qi qi qi+i ...] T ) such that a modified state-transition-matrix integrator can be used to march 
the coupled fluid-structural system forward in time. At the beginning of each time-step, the incremental 
structural deflections are calculated using the modal velocities and generalized aerodynamic forces. Using a 
deforming mesh, the mesh points are moved so that the inner mesh boundaries conform to the new deformed 
shape of the structure while the far-held boundaries are held stationary. Using a geometric conservation law, 
the fluid how through the updated mesh is obtained and used to calculate the generalized aerodynamic forces 
acting on the structure through the next time-step. Thus, a tight coupling of the how and the structure is 
implemented through the generalized aerodynamic forces. Thus, a time-history of the modal displacements, 
modal velocities and generalized forces is obtained. 

The aeroelastic capabilities of CFL3D, based on this modal response approach for obtaining the flutter 
boundary, have been partially validated for the transonic regime for the first AGARD standard aeroelastic 
configuration for dynamic response, Wing 445.6. The results of flutter calculations using Euler aerodynamics 
are given in Ref. 30 and those using Navier-Stokes aerodynamics are given in Ref. 31. However, these 
calculations were limited to the transonic regime. 


C. Computational Methods for Fluid-Structure Coupling 

Prediction of the dynamic response of a hexible structure in a fluid requires the simultaneous solutions of 
the equations of motion of the structure and the huid. In order to impose the kinematic boundary con- 
ditions on the fluid mesh at the new time step, the location and velocity of the fluid-structure boundary 
must first be known. This requires the solution of the entire system of equations for the structure, a task 
which cannot be carried out until the current surface pressure is known, which depends on the solution 
of the fluid domain and thus also on the unknown boundary conditions during the current time step. In 
addition, the discretized model of the structure uses a Lagrangian approach by following a point located on 
the structure over time, while the discretized model of the fluid uses an Eulerian approach by computing 
the flow quantities at a specific location in space over time. Therefore, accurate coupling of the two systems 
is a fairly complicated endeavor. A straightforward approach to the solution of the coupled fluid-structure 
system requires changing the fluid grid at each time-step, which is computationally very expensive. There- 
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fore, several different approaches have emerged as alternatives to partial regridding in transient aeroelastic 
computations, among them being dynamic meshes, 32 the space-time formulation, 33-35 the Arbitrary/Mixed 
Eulerian-Lagrangian formulation, 36,3 ' the multiple-field formulation, 38,39 the transpiration method, 19,40 the 
subgrid/TFI method, 41 and the modified spring analogy. 42 A detailed description of these methods is pro- 
vided in Ref. 43. However, since the CFL3D code used here implements fluid-structure coupling using the 
subgrid/TFI method and modified spring analogy, these methods are described in detail below. 

The subgrid/TFI method is an algebraic mesh deformation method in which surface movement is trans- 
mitted with an exponential decay into the mesh interior. The motion of selected slave points, chosen across a 
grid at constant index intervals, is tied to the motion of the nearest surface (or master) point. The exponen- 
tial decay function uses distance between the slave and nearest surface point so that motion of the surface 
is transmitted nearly undiminished to nearby slave points. On block faces, intervening mesh points are 
updated using transfinite interpolation (TFI). Once the intermediate mesh between slave points is updated 
on block faces, block interiors are updated using a volume TFI step. When smoothing of the deformed grid 
is necessary, several iterations of a modified spring analogy are employed. This scheme is a modification of 
the spring analogy 32 using axial spring stiffness. Spring stiffness in the mesh interior is controlled by the 
spacing of the appropriate boundary grid points. Note that the axial stiffness approach results in smoothing 
of the mesh and also allows adaptation based on the flow solution. Furthermore, the problem of grid collapse 
around convex surfaces is handled by selectively increasing/decreasing stiffness based on surface curvature. 
Finally, the results generated using CFL3D in this study employ an improved scheme that ensures the gen- 
erated control points at block boundaries are coincident and also checks block boundaries for separation 
during the mesh deformation at each time step. More details for this improved scheme will be available in 
Ref. 44. 

D. General Overview of the Solution Process 

The solution of the computational aeroelasticity problem used in the present study is shown in Fig. 2. First, 
the vehicle geometry is created using CAD software, and from this geometry a mesh generator is used to 
create a structured mesh for the flow domain around the body. In parallel, an unstructured mesh is generated 
that represents the surface of the finite element model of the vehicle. This structural model is further refined 
with internal stiffeners and mass elements to represent a complete vehicle, and is used to obtain the free 
vibration modes of the vehicle using MSC NASTRAN. Subsequently, the fluid mesh is used to compute the 
flow around the rigid vehicle using a CFD solver, which consists of the CFL3D code developed by NASA 
Langley Research Center. In order to generate the aeroelastic input for CFL3D, the modal deformation at 
each surface gridpoint in the fluid mesh is obtained by using cubic interpolation (in MATLAB) from the 
finite element structural model, for each structural mode. Using the flow solution as an initial condition, 
and the interpolated modal data as additional input, an aeroelastic equilibrium state is obtained for the 
flexible vehicle. For a geometry with vertical symmetry at zero angle of attack, such as a double-wedge 
airfoil, the equilibrium state is the same as the undeflected state. Next, the structure is perturbed in one 
or more of its modes by an initial modal velocity condition, and the transient response of the structure 
is obtained. To determine the flutter conditions at a given altitude, aeroelastic transients are computed at 
several Mach numbers and the corresponding dynamic pressures. The frequency and damping characteristics 
of the transient response for a given flight condition and vehicle configuration can be determined using some 
time domain damping and frequency identification method. This approach applied at the same altitude and 
vehicle configuration for a range of Mach numbers results in a series of damping values for the system. The 
flutter Mach number can be estimated from this series by interpolating the damping data points to identify 
the value of the Mach number at which the damping is zero. 

E. Time Domain Damping and Frequency Identification 

As mentioned earlier, the aeroelastic behavior of a configuration is obtained from CFL3D as a set of aeroelas- 
tic transient responses; from these the damping and frequency characteristics of the system are determined 
in the time domain. Two commonly used time domain identification methods are the moving block approach 
(MBA) 45,46 and the least squares curve fitting method (LSCFM). 4 ' A more recent method performs system 
identification of an Auto-Regressive Moving Average (ARMA) model of the system. 48-50 

The MBA, originally developed in the early 1970’s, uses the finite Fourier transform of a discretely 
sampled transient signal to identify modal damping and frequency. Consider a transient signal given by: 
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where 


q.i(t) = Aie CiUnit sin{uit + <j> a i) 


( 5 ) 


a;? = u&(l-C?) 

Then, the finite Fourier transform of this signal is: 

[Tt+Tr 

Q'i(ui,T t )= / 4,e" Ci “" i ‘sm(w l ( + 

J Tt 


i)e~ iUit dt 


( 6 ) 

( 7 ) 


where Qi( w »i T t) * s a function of r t at the analysis frequency, u>i. The moving block function is defined as the 
amplitude of Eq. (7). When the damping is small, i.e. Q <C 1, the natural logarithm of the moving block 
function is given by: 


InQtfai, T t ) = -Ci^iTt + - (i sin{2(u>iT t + 4> a i)} + constant ( 8 ) 

Note that Eq. (8) has a linear component with a slope of and an oscillatory component at twice the 

analysis frequency, uj % . 

The method is applied by using a Fast Fourier Transform (FFT) algorithm on the entire transient response 
to calculate the frequency spectrum of the signal, from which a frequency of interest, or analysis frequency, 
is selected. Next, Eq. (8) is calculated on a block of the data, where N b < N, at r t = 0. This is repeated 
N — N b times, after shifting the block one sample point at a time, i.e. r t = nAt for n = 0, 1, 2, ..., N — Nt,. 
The slope of Eq. (8) is calculated from the resulting curve using a least-squares curve fit, and the damping 
is calculated by dividing this slope by the analysis frequency. 

A critical step with this method is the selection of an analysis frequency from the frequency spectrum 
of a given response. Since the aeroelastic response of the system is produced as a set of transient modal 
amplitudes, there are two options. One option is to calculate the frequency spectrum of the response from the 
transient displacement, Eq. (1), of a point on the structure. The analysis frequencies of the displacement will 
correspond to local maxima in the frequency spectrum. Image processing, therefore, is required, either by 
the user, or by an image processing subroutine, to locate the frequencies. The second option for locating the 
frequencies is to calculate the frequency spectrum of each transient modal amplitude individually. However, 
since the modal amplitudes are coupled through the generalized forces, the frequency spectrum of each 
modal amplitude will contain multiple frequencies. The analysis frequency, therefore, should be selected as 
the global maximum from the frequency spectrum of each modal amplitude. Note, however, that as flutter 
is approached, the frequencies of some modal amplitudes can become dominant throughout the system. 
When this occurs, selecting the global maximum from the frequency spectrum results in sudden shifts of the 
damping and frequency of some of the transients, since the MBA is tracking a different component. While 
information of some of the modes is lost when this occurs, this option provides the analyst insight as to how 
the modes interact as the flutter Mach number is approached. 

There are a two main issues with the MBA for damping and frequency identification. In cases where 
a particular mode damps out rapidly, if the response record is too long, the moving block function will no 
longer form a straight line, and therefore the damping, which is calculated from a linear least squares fit of 
the MBA output data, will be incorrect. Since the amount of damping in each mode is not known before a 
response is processed, this issue introduces uncertainty in appropriate record lengths for accurate damping 
estimates. A second issue is the inability to identify peaks in the frequency spectrum when there are closely 
spaced modal frequencies. 

The LSCFM, 4 ' also developed in the 1970’s, performs a least squares curve fit of a transient response to 
exponential functions in order to identify modal damping and frequency. 

Consider the following complex exponential function: 

rim 

f(t ) = ao + e -Ci^nit { ai cos(uiit) + bisin(uiit)} (9) 

i= 1 

where n m represents the number of modes in the signal. The modal damping and frequency of a transient 
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response can be calculated by minimizing the squared-error difference between the output fit f(t k ) and the 
input transient response w(t k ), where tk is the k th discrete time sample of the response. This error is given 
by: 

N 

E = - W (t k )} 2 (10) 

fc= 1 

Note that this is a nonlinear least squares problem. One method of solution is to provide initial “guess” 
values for Q and Cu\, and then solve the linear least squares problem for ag, a,i, and © The curve fit is then 
calculated by searching for the Q and u>i which minimize Eq. (10). At each search step the linear coefficients 
(ao, a-i, and bi) are updated, and the error is recomputed, using the current Q and w,. The search for Q and 
u>i in this study is performed using the subroutine FMINSEARCH in MATLAB©. 

The weakness of this method is the computational time required to find the a o, a,;, and bi, which minimize 
Eq. (10), particularly when more than two modes are present in the aeroelastic system. 

The ARMA system identification method for damping and frequency identification 48-50 is based on a 
single-input single-output deterministic ARMA model of the aeroelastic system, with 2 n m Auto-Regressive 
(AR) and one Moving Average (MA) coefficients. The model has the following form: 


2n rt 


Qk + O’iQk-i — biSk-1 


( 11 ) 


where 2 n m AR coefficients are used to determine the aeroelastic system damping and frequencies, and one 
MA coefficient is sufficient for identifying an aeroelastic static offset. To identify the damping and frequency, 
let the transient aeroelastic response be represented by the following AR model: 


2 n„ 


Qk + ^2 a i < ik—i = 0 


i= 1 


This can be written in state-space form as: 


( 12 ) 
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The state vector, { X p } k , is defined as 


(14) 
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with 



qk = 

— ai9fc-i + h\{k — 1) 


h\(k) = 

Q j 2Qk— 1 + h 2 {k — 1) 

^2n n 

,-2 (fc) = 

—a-2n m -iqk-l + ©n„ 

^2n„ 

,-iO) = 

— d2n m qk- 1 


( 17 ) 


Note that the state-space description of the AR model in Eq. (13) is in observer form and is completely 
observable. The aeroelastic system damping and frequencies are determined from the eigenvalues of the 
estimated matrix, {A p }, given in Eq. (14). These can be written as: 


Aj = r, ■ i Sj 
Aj+n m = r j ~ i s j 


(18) 


where j = 1, 2, • ■ ■ , n m . Note that the AR model is in the discrete time domain, therefore, the aeroelastic 
modal damping and frequencies in the continuous time domain are given by: 


0 = sk l °9e{r 2 j + s j) 

ujj = 7^r tan ~ 1 ~ 


(19) 


In this study, the AR coefficients are determined from the aeroelastic responses using the subroutine AR 
from the system identification toolbox in MATLAB©. Note that the sample times, used to calculate the 
AR coefficients from the response data, were set to: 


T e = 


2(w n ) 


n J max 


( 20 ) 


where (uj n ) max is the maximum system natural frequency in radians/second. 


F. Computational Model for the Three-Dimensional Low Aspect Ratio Wing 

Few studies 51,52 have been carried out that validate CFL3D for the hypersonic regime. In part, this is 
because CFL3D was not originally designed for modeling hypersonic flow, and does not contain the necessary 
ingredients for computing real-gas effects. However, depending on the cases considered, and the most relevant 
flow variables, sufficiently accurate solutions can be obtained while ignoring real gas effects. References 51 
and 52 compared the pressure distributions and static and dynamic stability derivatives of cones and ogive- 
cylinder bodies obtained using CFL3D with results obtained using a unified hypersonic/supersonic panel 
method in the range 3 < M < 6. However, no attempt was made in these studies to compare aeroelastic 
stability boundaries obtained with CFL3D. 

A systematic validation of the CFL3D code for flutter analysis in the hypersonic regime has never 
been undertaken. To validate aeroelastic stability boundaries computed with CFL3D, it is important to 
identify and use simple configurations for which aeroelastic stability boundaries can be computed using 
alternative approaches. Specifically, generating results for simple three-dimensional configurations using 
Euler and Navier-Stokes unsteady aerodynamic loads, and comparing them with results obtained using an 
independently developed aeroelastic code based on piston theory, represents a validation of CFL3D in the 
hypersonic regime. 

The computational model for the low aspect ratio wing, shown in Figs. 3 and 4, was developed from a 
study 26 of grid configurations, where the appropriate computational domain and mesh resolution required 
for the hypersonic aeroelastic analysis of the low aspect ratio wing were determined. In Ref. 26, cell 
distribution and grid resolution were varied to construct a grid for the hypersonic flow regime that accurately 
and efficiently captured the flow. Each span section plane of the computational domain extends one-half 
chord-lengths downstream. The boundary of the grid surrounding the wing, from the leading edge to mid- 
chord, extends to a distance 10% beyond the shock that forms at M = 5.0. The computational domain in 
the spanwise direction also extends beyond the tip of the wing by 35% of the semi-span length. Furthermore, 
the grid is tapered, in all three dimensions, so as to be compatible with the geometric taper of the wing. 
With 0.63x10® cells, this grid is a 57x353x33 C-gricl with 353 points around the wing and its wake (289 
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points on the wing surface), 57 points extending spanwise from the root (49 points on the wing on the wing 
surface), and 33 points extending radially outward from the surface. 

G. Computational Model for the Generic Hypersonic Vehicle 

In order to carry out aeroelastic simulations of the generic vehicle, shown in Fig. 5, using Euler aerodynamics, 
various meshes were constructed for the computational domain shown in Fig. 6. Developing a suitable mesh 
for use with CFL3D requires balancing between various conflicting constraints. One of the main factors 
governing the usefulness of the mesh is the overall size of the mesh, which directly affects runtime. Having 
fixed the overall mesh size, obtaining accurate results requires having optimum placements for the internal 
gridpoints. This is done by having the maximum possible number of gridpoints lie within the shocks present 
around the vehicle. The midplane of the selected mesh is shown in Fig. 7, while Fig. 8 shows the gridpoints 
around the typical cross-section of the generic vehicle. On the surface of the generic vehicle, the gridpoints 
are again finely spaced close to the surface of the vehicle. 

By varying the nondimensional grid-spacing, representing the distance of the closest grid-point in the 
normal direction, different pressure distributions on the surface of the vehicle’s canted fin, at 75% of span, 
were obtained in order to assess convergence. The grid size of 3 x 10® gridpoints, the gridpoint placement 
on the surface and all other grid parameters were kept constant for these simulations. As seen in Table 
1, the error in the moment acting on the typical section of the canted fin was reduced with decreasing 
Afc. Reference 26 indicates that the maximum error of 0.28% is small enough to expect converged results. 
However, aeroelastic simulations could not be carried out for less than 5 x 10 -7 , as this produces mesh 
entanglement. When this occurs, cell boundaries penetrate through adjacent cell volumes. In CFL3D, this 
generates negative cell volumes in the region surrounding the entanglement. Reducing the time-step enables 
the use of slightly smaller grid spacing, but with a significant increase in computational time. Hence, the 
time-step size was selected to be At = 0.35 x 10 -3 sec, which allows accurate calculation of the aerodynamic 
forces, provides at least 50 time-steps per cycle of the highest frequency mode, and avoids mesh entanglement 
for this value of A^. The above mesh was referred to as the “fine” mesh, and for comparison, an additional 
mesh having 4 x 10 5 gridpoints was also created. This mesh is referred to as the “coarse” mesh, and has a 
nondimensional grid-spacing at the surface in the coarse grid of A& = 6.5 x 1(R 6 . 

The coarse mesh is a 65x177x33 C-grid with 177 points around the vehicle and its wake, 65 points 
spanning the grid from side to side, and 33 points extending radially outward from the vehicle surface. 
Of these, 145 points were wrapped around the fuselage in the streamwise direction, while 33 points were 
distributed on the vehicle fuselage from side to side. The fin surfaces were discretized using 9 grid points 
in the span-wise direction and 17 grid points from leading edge to trailing edge. The refined mesh is a 
129x355x65 C-grid with 355 points around the vehicle and its wake, 129 points spanning the grid from 
side to side, and 129 points extending radially outward from the vehicle surface. Of these, 289 points were 
wrapped around the fuselage in the streamwise direction, while 65 points were distributed on the vehicle 
fuselage from side to side. The fin surfaces were discretized using 17 grid points in the span-wise direction 
and 33 grid points from leading edge to trailing edge. In both cases, the computational domain was extended 
0.5 vehicle-length upstream and 1.5 chord lengths downstream. On either side, the domain extended 0.5 
vehicle-lengths spanwise. The distance to the upper and lower boundaries was selected to be slightly larger 
than the shocks forming at M= 5.0, with the shocks from the nose of the vehicle flowing out of the domain 
boundaries without reflection. 


H. Piston Theory Aerodynamics 

Piston theory is a simple inviscid unsteady aerodynamic theory, that has been used extensively in supersonic 
and hypersonic aeroelasticity. It provides a point-function relationship between the local pressure on the 
surface of the vehicle and the component of fluid velocity normal to the moving surface. 53,54 The derivation 
utilizes the isentropic “simple wave” expression for the pressure on the surface of a moving piston, 


where 


P{x,t) = f 1 + 7 - 1 Vn \ 

Poo \ ^ ^oo ) 


V n 


dZ(x,y,t) dZ(x,y,t) 

dt dx 


( 21 ) 


( 22 ) 
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The expression for piston theory is based on a binomial expansion of Eq. (21), where the order of 

V 

the expansion is determined by the ratio of — Reference 54 suggested a third-order expansion, since it 

produced the smallest error of the various orders of expansion used when compared to the limiting values 
of pressure, namely the “simple wave” and “shock expansion” solutions. The third-order expansion of Eq. 
( 21 ) yields third order piston theory: 


p(x, t)-p oa = Poo 


7 - 


«n ,7(7 +1) / 7 ( 7 +!) ( V, X 3 


4 V a oo/ 12 

Combining Eqns. 1, 4, 22, and 23 yields the generalized forces obtained from piston theory. 


(23) 


I. Refined Aerothermoelastic Model 

A realistic aeroelastic model for the hypersonic regime must include aerodynamic heating effects. Aerody- 
namic heating significantly alters the flow properties , 55 degrades the material properties and also introduces 
thermal stresses . 56-58 Aerodynamic heating of the flow surrounding the vehicle leads to significantly differ- 
ent thermodynamic and transport properties, high heat-transfer rates, variable 7 , possible ionization, and 
nonadiabatic effects from radiation . 55,56 Thermal stresses can arise from rapidly changing conditions of heat 
input where time lags are involved, or from equilibrium conditions of non-uniform temperature distribu- 
tion . 57,58 Commonly, the heated structure has lowered stiffness due to material degradation and thermal 
stresses, which manifests themselves as a reduction in frequencies . 57-59 

Previously, the authors of this study performed an exploratory aerothermoelastic analysis on a low aspect 
ratio cantilevered wing . 26 As a first approximation to the aerothermoelastic system, the heating distribution 
of the system was assumed to be a parabolic function of the chord. Furthermore, steady-state conditions 
were assumed, resulting in a unique temperature distribution for each Mach number. An exact treatment of 
aerothermoelasticity requires the coupling of the unsteady heat transfer problem with the aeroelastic problem 
based on a Navier-Stokes solution of the unsteady airloads, which results in time dependent temperature 
distributions. This implies time dependent free vibration characteristics of a structure, at a given Mach 
number, as it is heated. 

The heat transfer between the fluid and the structure, schematically depicted in Fig. 9, is determined 
from an energy balance of the heat fluxes at the wall of the structure : 60,61 


where 


Qaero — Qrad T Qcond Qstrd 


Qaero — ^ht i^-AW ) 


Qrad — CFcT] 


W 


dT\ 

Qcond — ( ^7^ J 

dT 1 


Qstrd PwCpw^w 


W 


dt 


(24) 

(25) 

(26) 

(27) 

(28) 


The heat transfer represents a balance at the wall between the convective heating by the fluid ( q aer o ) 
and heat loss due to conduction into the structure ( q CO nd ), radiation out to space ( q ra d ), and energy stored 
in the wall ( q s trd )■ The heat transfer problem is driven by Eq. 25, or more specifically the adiabatic wall 
temperature, T A yy, which is a function of the surface geometry and free stream conditions. Note that at 
steady-state, 

Qaero — Qrad (29) 

Therefore, using Eqs. 25 and 26, 

— h ht T R — h ht T AW = 0 (30) 

where, given hht, T A Wi cr, and e, the radiation equilibrium temperature at the wall, Tr, can be determined 
by solving a 4th order algebraic equation. 
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The solution procedure for the refined aerothermoelastic system is depicted in Fig. 10. Comparing 
Figs. 2 and 10, the modification introduced by the aerothermoelastic solution is due to the fact that 
the aerodynamic heating information, namely Taw and h^t from Eq. 25, obtained on rigid body CFD 
computations, is introduced into the finite element analysis of the system. Specifically, Taw is calculated by 
specifying an adiabatic wall as a boundary condition for the CFD solution to the Navier-Stokes equations, 
and hht is calculated from a separate Navier-Stokes solution where the wall temperature boundary condition 
is set to free stream temperature. Setting a free stream, or “cold”, wall boundary condition to calculate h^t 
provides a conservative estimate of the aerodynamic heating. The transient temperature distribution in the 
structure is determined from a heat transfer analysis using finite element analysis, and the free vibration 
frequency and mode shapes of the transiently heated structure are then calculated at each desired point in 
time. This process is depicted in Fig. 11. 

Note that this is not an exact treatment of the aerothermoelastic system, since the temperature dis- 
tribution is not computed at each time step of the aeroelastic calculation procedure. It is, however, a 
reasonable approximation since the heat loads vary slowly with time compared to the generalized forces and 
motion of the system. Furthermore, it is computationally prohibitive to generate an aeroelastic computation 
continuously for the time scales involved to heat the structure. 

The aerothermodynamic quantities in this study were determined using a 2.1x10® cell grid using CFL3D 
to solve the Navier-Stokes equations. The heat transfer analysis was performed using the “Transient Analysis 
Solution” (Sol 159) 62 in MSC.NASTRAN with the “HEAT” option selected. The effect of thermal radiation 
was included. The heated modes and frequencies of the wing were determined using the “Nonlinear Statics 
Solution” in MSC.NASTRAN (Sol 106) 62 with the “NORMAL MODES” option selected. The analysis 
includes both the effect of material property degradations and thermal stresses. 

III. Results and Discussion 

Comparison of aeroelastic stability results based on Euler, Navier-Stokes and piston theory in Refs. 
24-26 provided useful guidelines regarding the importance of viscosity and aerodynamic heating, and the 
effectiveness of piston theory in approximating the aeroelastic behavior of simple two- and three-dimensional 
lifting surfaces. A similar approach is pursued here for three-dimensional low aspect ratio wings. First, 
however, the effectiveness of three time domain damping and frequency identification techniques is studied 
using aeroelastic data generated from a simple typical section analysis using unsteady Navier-Stokes aerody- 
namics. 25 Next, the low aspect ratio wing is used to compare the aeroelastic behavior generated using third 
order piston theory, Euler and Navier-Stokes aerodynamics. Furthermore, the sensitivity of the aeroelastic 
behavior to several parameters related to time stepping is studied using Euler aerodynamics. Next, the 
aerothermoelastic behavior of a low aspect ratio wing is studied using the refined aerothermoelastic model. 
Finally, the aeroelastic behavior of a generic hypersonic vehicle is examined using third order piston theory 
and Euler aerodynamics at various altitudes. This configuration is also used to study the sensitivity of the 
aeroelastic behavior to grid resolution using Euler aerodynamics. 

A. Comparison of Time Domain Frequency and Damping Identification Techniques 

The three identification methods (MBA, LSCFM, ARMA) have been previously compared, in Ref. 49, using 
several damped free oscillation examples, as well as subsonic and transonic aeroelastic data. The ARMA 
method was found to be superior for the cases considered. 49 Specifically, the ARMA model calculated the 
damping and frequency of modes with only 0.01% separation in frequencies. Furthermore, the ARMA model 
required substantially smaller time records than either the MBA or LSCFM to accurately extract system 
damping and frequencies. For computational aeroelasticity, this is an important consideration when choosing 
a time domain flutter identification technique since the cost of a computational aeroelastic simulation is 
largely determined by the size of the response record required to calculate system damping. 

To compare these three methods, each was used to calculate the frequency and damping for a simple 
typical section aeroelastic system, shown in Fig. 12. The parameters describing this configuration are listed 
in Table 2. The aeroelastic responses were generated by CFL3D using unsteady Navier-Stokes aerodynamics 
at an altitude of 40,000 ft. The computational domain used for the CFD computations of the typical section 
was created from the 75% span section of the low aspect ratio wing computational domain. The aeroelastic 
transients were calculated for 4000 time steps, corresponding to 1 second of response. A comparison of each 
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method is given in Fig. 13. For the results labeled ARMA, MBA, and LSCFM, the damping and frequency 
of the system were calculated from the pitch response of the system, where: 


h(t) 

a(t) 



(31) 


The results labeled MBA2 were generated using the individual aeroelastic transients, i.e. {< 7 i(i),<? 2 (t)}- 
The MBA results, therefore, were generated using the MATLAB© subroutine IMREGIAONALMAX to 
locate the analysis frequencies from the frequency spectrum of the pitch response. The MBA2 results were 
generated using the global maximum from the frequency spectrum of each individual transient as the analysis 
frequencies. 

It is evident that all three methods calculate similar damping and frequency estimates for the typical 
section responses. In particular, each method predicts the same flutter Mach number, Mf = 11.9, when 
the damping becomes zero. Furthermore, the ARMA, LSCFM, and MB A/MB A2 methods predict the same 
aeroelastic behavior, over the range of Mach numbers considered, for the second mode. The behavior of the 
first mode, i.e. lower frequency mode, is calculated equivalently using both ARMA and LSCFM. The MBA 
and MBA2 methods show the behavior of the first mode as approaching the second mode as the flutter 
Mach number is approached. Specifically, for the MBA2 method, the aeroelastic behavior of the first mode 
is similar to the second mode from M = 11.0 up to flutter. Frequency spectrum plots of the transients, 
shown in Figs. 14 and 15, reveal that the dominant frequency (i.e. global maximum of FFT) of the first 
modal transient shifts from the first modal frequency to the second modal frequency between M = 10.0 
and M = 11.0. This indicates that the second modal frequency dominates both the first and second modal 
transients as flutter is approached. For the MBA method, the damping of the first mode approaches the 
damping of the second mode from M = 11.0 to M = 11.75. Beyond M = 11.75, no analysis frequency is 
available for the first mode. Figures 16 and 17 illustrate that at M = 11.0, two local peaks are present in 
the frequency spectrum of the pitch response. However, at M = 11.75, the modal frequencies have coalesced 
to a point where the local peak corresponding to the first modal frequency is no longer identifiable by the 
image processing routine, IMREGIONALMAX. This illustrates the failure of this method to identify closely 
spaced analysis frequencies. 

The performance of each method as the time record of the response is reduced is illustrated in Figs. 
18 - 21. It is evident from these results that, for this aeroelastic system, both the ARMA method and 
the LSCFM are capable of computing the aeroelastic behavior accurately when the time record is reduced 
substantially. In this case, the aeroelastic behavior is calculated from these two methods with only 250 
time steps, corresponding to 0.06 seconds. However, the performance of the MBA and MBA2 methods 
diminishes as the time record is shortened. These methods require at least 0.25 seconds of response data to 
provide reasonable estimations for the damping and frequency. Furthermore, the estimated behavior, using 
ARMA and LSCFM, is invariant when the response record is shortened, while there is some variation in the 
MBA/MBA2 results as the record length is changed. In both the MBA/MBA2 results, the flutter Mach 
number is slightly reduced, and the behavior before flutter varies, as the record length is shortened. 

In the cases considered here, the ARMA method is superior to both the MBA and LSCFM. It produces 
damping and frequency estimates from time records that are 75% shorter than those required for the MBA. 
Furthermore, it is superior to the LSCFM in terms of computational efficiency. The ARMA method required 
~3 seconds on a 3.06 GHz Xeon processor, while, depending on the size of the time record used, the LSCFM 
required up to 5 minutes of CPU time. Furthermore, the time required for the ARMA method to generate 
damping and frequency estimations is relatively independent on the number of modes tracked in the system, 
while the time required for the LSCFM increases with the number of modes tracked. 


B. Aeroelastic Behavior of a Three-Dimensional Low Aspect Ratio Wing 

The parameters describing the physical properties of the three-dimensional low aspect ratio wing are provided 
in Table 3 and Fig. 22. The natural frequencies and modes, shown in Fig. 23, were determined by comparing 
them with the bending and torsional frequencies and total mass of a wing that resembles the Lockheed F-104 
wing. For the Euler and Navier-Stokes computations, the time step size was set to accommodate 50 steps 
per cycle of the highest frequency mode, which corresponds to At = 0.25xl0~ 3 seconds. This resulted in a 
smaller time step than that required for modeling the unsteady aerodynamic loads. 25 
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The sensitivity of the aeroelastic behavior of the wing to the number of modes is illustrated in Table 4. It 
is evident that increasing the number of modes from 5 to 8 has little effect on the aeroelastic behavior when 
third order piston theory aerodynamics is used. However, it is also clear that at least five modes are needed 
to capture the aeroelastic behavior. The aeroelastic behavior of the low aspect ratio wing using third order 
piston theory with 5 modes, at 40,000 ft, is shown in Fig. 24, where Mf = 13.4. Using the frequency of the 
flutter mode, this corresponds to = 0.14, while for Mode 5, k u = 0.51. For this case, as the Mach number 
is increased, there is a coalescence in frequency for the first bending (Mode 1) and first torsion (Mode 2) 
modes, while the damping of the first mode approaches zero. Also, as the damping of Mode 1 decreases, 
the damping of Mode 2 increases. The damping and frequency of the Modes 3-5, are relatively constant 
with increasing Mach number. Note that the ARMA system identification method was used to calculate the 
damping and frequency of the system. 

In order to solve the unsteady Euler and Navier-Stokes equations, several temporal parameters must be 
set in CFL3D so as to conduct efficient time-accurate calculations. Specifically, CFL3D employs a “pseudo” 
time step 27 in order to reduce linearization and factorization errors in the time-accurate computations. The 
“pseudo” time step is implemented by specifying a number of sub-iterations, as well as the “pseudo” time step 
size, which is governed by CFL T . Generally speaking, decreasing the global time step, reduces the number 
of sub-iterations required to achieve an accurate result, while increasing the grid refinement increases the 
number of sub-iterations required. 27 The effect of sub-iterations and CFL T on the hypersonic aeroelastic 
behavior of the low aspect ratio wing is illustrated in Fig. 25 using a 0.63x10 s cell grid. These results were 
generated at 40,000 ft, for M= 12.0, using Euler aerodynamics. This Mach number and altitude was chosen 
since it was expected, based on the piston theory results, to result in aeroelastic behavior that was relatively 
close to, but not above, the flutter boundary. It is evident from Fig. 25 that the number of sub-iterations 
used has an effect on the aeroelastic behavior of the wing. While the frequencies of all the modes remain 
constant, the damping varies with increasing sub-iterations. In particular, the first and second modes are 
substantially affected by the number of sub-iterations. The most significant impact is the relative magnitude 
of the damping, where the first and second modes cross each other depending on the number of sub-iterations 
used. For CFL T = 1.0, the damping of the second mode is greater than the damping of the first mode, based 
on 25 - 45 sub-iterations. Furthermore, the damping of all the modes remains almost constant in this range 
of sub- iterations. Also, when CFL T = 5.0, which represents a larger step in “pseudo-time” per sub-iteration, 
the same behavior is observed, but for a fewer number of sub-iterations. In the CFL T = 5.0 case, when the 
number of sub-iterations varies between 20 - 35, the behavior is relatively unchanged. Based on these results, 
CFL t = 5.0 with 20 sub-iterations is selected since it provides the most efficient and accurate results when 
Euler aerodynamics are used. As illustrated in Fig. 26, the flutter Mach number of the low aspect ratio 
wing at 40,000 ft is Mf = 13.7 (0.15 < k u < 0.49) when Euler aerodynamics are used. Note that this is 
only 2% higher than the flutter Mach number predicted using third order piston theory. Also, similar to the 
piston theory results, there is a coalescence in frequency between the first bending (Mode 1) and first torsion 
(Mode 2) modes, while the damping of the first mode approaches zero. Furthermore, the same divergence 
of damping in the first two modes is observed as the flutter point is approached. Again, similar to the 
piston theory results, the damping and frequency of Modes 3-5 remains nearly constant with increasing 
Mach number. In order to asses the effect of grid resolution on the flutter boundary, these computations 
were repeated on a 0.27x10 s cell grid. These results are shown in Fig. 27, where Mf = 13.8 (0.15 < k u < 
0.51). It is evident that coarsening the mesh by almost 60% does not significantly alter the flutter boundary 
at this altitude. However, it is important to note that as the flutter Mach number is increased, such as 
due to an increase in altitude, the shock will lie closer to the surface of the wing. Therefore, as the shock 
approaches the surface, grid resolution may become significant since the number of the points in the shock 
layer decreases. Note that the results generated using the 0.63x10 s cell grid required 6.5 hours of CPU time 
using 12 AMD Athlon 2000MP CPUs, where 500 time steps were used to generate 0.125 seconds of transient 
response. The ARMA system identification method was used to calculate the damping and frequency of the 
system. 

Using the number of sub-iterations required for the Euler results as a guideline, the aeroelastic behavior 
of the low aspect ratio wing was determined for Navier-Stokes aerodynamics. For an initial estimate of 
the aeroelastic behavior using Navier-Stokes aerodynamics, CFL T = 1.0 with 50 sub-iterations was chosen. 
The results, shown in Fig. 28, are similar to those generated using third order piston theory and Euler 
aerodynamics. In this case, Mf = 13.6 (0.14 < k u < 0.48), which is only 1% lower than that predicted 
using Euler aerodynamics, and is 1.5% higher than that predicted using third order piston theory. As with 
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the piston theory and Euler results, there is a coalescence in frequency, and a divergence in damping, of the 
first two modes, as the flutter Mach number is approached. Furthermore, the behavior of the higher modes 
is similar, where there is little change in the damping and frequency as the Mach number is increased. A 
notable difference between these results and the inviscicl results is the gradual approach of the first mode 
to zero damping. Both the piston theory and Euler results had a more sudden change in damping at the 
flutter Mach number. In order to asses the effect of CFL T on the flutter boundary, these computations were 
repeated using CFL T = 5.0. These results are shown in Fig. 29, where M/= 13.65 (0.14 < k u < 0.50). 
The small change in the aeroelastic behavior and flutter Mach number due to the substantial increase in 
CFL t suggests that these results are converged to the correct solution. However, these aeroelastic results 
generated using Navier-Stokes aerodynamics should still be considered preliminary since the effect of sub- 
iterations has not been carefully studied. Also, it has been shown 26 that grids on the order of 2xl0 6 cells are 
required to converge the drag coefficient in steady-state computations. The effect of grid resolution on the 
hypersonic aeroelastic behavior of the wing is not known, and therefore needs further attention. These results 
appear reasonable, however, since it has been shown in Ref. 25 there are only small differences between the 
aeroelastic behavior of double-wedge typical section when either Euler or Navier-Stokes aerodynamics are 
used. Note that these results were generated on a 0.63xl0 6 cell grid, requiring 16.5 hours of CPU time 
using 12 AMD Athlon 2000MP CPUs, where 500 time steps were used to generate 0.125 seconds of transient 
response. The ARMA system identification method was used to calculate damping and frequency of the 
system. 

Note that the reduced frequencies, k^, at flutter are small irregardless of aerodynamic method used. 
This explains partially the good agreement between the Euler/Navier-Stokes results and the piston theory 
results. It is interesting to note that the aeroelastic behavior predicted using piston theory is conservative 
compared to the aeroelastic behavior predicted using Euler and Navier-Stokes aerodynamics. This is likely 
due to both three-dimensional flow effects which are not accounted for by piston theory, since it is strictly a 
point function relationship between the body and the pressure, and also, the small differences between the 
viscous and inviscid results for this configuration. 

C. Aerothermoelastic Behavior of a Three-Dimensional Low Aspect Ratio Wing 

The low aspect ratio cantilevered wing is a convenient example for studying the effects of aerodynamic 
heating on an aeroelastic system. For such a configuration the restrained warping at the root of the wing 
will induce thermal stresses which in turn will affect the torsional stiffness of the wing and modify its 
frequencies and mode shapes. The effect of warping restraint increases as the aspect ratio of a structure 
diminishes. 63 Early studies of this effect 64 indicate that by modeling a low aspect ratio wing as a plate, the 
effect of warping restraint is inherently included. More recently, this effect has been studied in the context 
of composites. 65,66 In these studies, it has been pointed out that warping restraint is not only important 
for low-aspect ratio metallic structures, but also for composite structures where the material proprieties are 
non-isotropic. Furthermore, it was shown that the warping stiffness of a cantilever plate is a function of both 
aspect ratio, and the ratio of bending/torsion stiffness. 1,65 This observation is important for hypersonic 
vehicles where structural properties are altered by aerodynamic heating. 

In Ref. 26, it was determined that the low-aspect ratio wing was prone to thermal buckling due to 
local deformations along the leading and trailing edges. Therefore, the wing was modified for this study, 
in order to reduce the susceptibility to thermal buckling. Specifically, the leading and trailing edges have 
been stiffened, and several vertical spars have been added throughout the entire wing. This modification 
introduced minor changes in the structural mass and first two free vibration frequencies. The mass and 
vibration frequencies of the original and modified low aspect ratio wings are compared in Table 5. The 
mode shapes of the modified wing are shown in Fig. 30. Note that only the fourth mode shape has changed 
significantly. For the original wing the fourth mode is an antisymmetric bending mode, and for the modified 
wing it is a symmetric bending mode. 

Since the low aspect ratio wing is constructed of aluminum, it must be augmented with a thermal protec- 
tion system (TPS) in order to provide reasonable heat transfer results in the high temperature environment 
of hypersonic flow. The TPS selected is based on Ref. 67. It consists of a 0.45mm RENE’ 41 metal heat 
shield, and a 3.8mm flexible MIN-K thermal insulation blanket between the wing and heat shield. The 
RENE’ 41 heat shield can withstand temperatures up to approximately 1500 K, and is assumed to have an 
emissivity, e = 0.85. 

The normalized radiation equilibrium wall temperature for several Mach numbers and altitudes are shown 

15 of 47 


American Institute of Aeronautics and Astronautics 


in Fig. 31. The radiation equilibrium wall temperatures, Tr, were normalized by the peak equilibrium 
temperature, (TR) maa; , which is listed below each case. These plots provide insight into the distribution 
in heating, as well as the magnitude of temperatures involved. In general, the heating is maximum at 
the leading edge, and minimum at the trailing edge. Also, note that as the Mach number is increased, 
the temperature distribution variation with chord increases. Furthermore, increasing the altitude, increases 
the amount of chordwise variation in temperature. These plot illustrate, that even for a wing augmented 
with a radiative heat shield, the temperatures at equilibrium in hypersonic flow will surpass the melting 
temperature of aluminum (^900 K) without active cooling. Furthermore, the temperatures present at M = 
10.0, at 80,000 ft and lower, will surpass the maximum temperature sustainable by the heat shield. 

Figure 32 illustrates the effect of transient heating on the free vibration frequencies of the modified wing 
at 100,0000 ft for M = 10.0, 15.0, and 20.0, at zero degree angle of attack. In each case, there is a reduction 
of the second, third, and fifth modal frequencies as the wing is heated. Also, there is a slight increase in the 
first frequency, and a more dramatic increase in the fourth modal frequency, which switches order with the 
fifth modal frequency. Also, in each case the wing buckles after a specific duration, namely 30 minutes, 18 
minutes, and 13 minutes for M = 10.0, M = 15.0, and M = 20.0 respectively. An interesting observation 
is the qualitative similarity between these three results. Each Mach number produces similar changes in 
frequency, where the largest difference in the results is the amount of time required for the changes to 
occur. This is likely due to increases in the heating rate due to increasing Mach number, while the heating 
distribution is relatively unchanged. In Fig. 33, the same results are presented, however, the frequencies 
are plotted as a function of a reference temperature, TR e j. The reference temperature was chosen as the 
leading edge temperature of the 75% span location. Note that there is only minimal differences in the three 
cases when the results are plotted as a function of temperature. This plot illustrates, that at 100,000 ft, the 
wing is susceptible to buckling, for 660 K < Tn e f < 680 K, irregardless of the Mach number. Furthermore, 
the free vibration characteristics of the heated structure at any Tfj e y, can be estimated from one of the 
cases, implying that an analyst can reduce computation time required to calculate the effect of aerodynamic 
heating on the structure. 

In practice, control surfaces on hypersonic vehicles will have a small angle of attack. It has been shown, 25 
that for the unheated case, a small angle of attack will not significantly alter the flutter boundary of a 
two dimensional typical section. When aerodynamic heating is considered, however, angle of attack is an 
important consideration since it introduces asymmetry in the temperature distribution between the upper 
and lower surfaces of a wing, and therefore introduces thermal stresses not present in zero degree angle of 
attack cases. The effect of angle of attack on the heated free vibration frequencies of the modified wing is 
shown in Fig. 34 at 100,000 ft for a s = 2°, and M = 20. It is clear that, in the heated case, a small of angle 
of attack has a significant impact on the free vibration frequencies of the wing. Specifically, there is not a 
smooth change in frequency as the wing is heated, nor is there a clear distinction of where or if buckling 
occurs. Furthermore, the frequencies are much lower in the angle of attack case for 600 K < Tn e f < 665 K. 

To gain insight into the effect of aerodynamic heating on the flutter boundary, the heated mode shapes 
of the wing were used to calculate the flutter boundary at 100,000 ft using third order piston theory aerody- 
namics. The heated configuration at a s = 2°, TR e / = 665 K was chosen since it appears from Fig. 34 that 
the structure is significantly weakened at this point due to heating, yet is not buckled. The flutter boundary 
of the heated wing and cold wing are compared in Fig. 35. The ascent trajectory of the X-43, which is 
similar to the ascent trajectory of the now defunct NASP vehicle, is also included in order to illustrate the 
proximity of the wing flutter boundary to typical operating conditions. Note that due to the extremely high 
flutter boundary of the wing at the higher altitudes, the flutter boundary is also calculated using unsteady 
Newtonian impact aerodynamics, 16,68 and used for comparison. At 100,000 ft there is a relatively small 
reduction in Mf of the cold wing ( Mf = 71.7) when a small amount of angle of attack ( Mf = 62.6) is 
introduced. Furthermore, at 100,000 ft, Mf is large implying that the piston theory results are inaccurate. 
However, they are conservative when compared to the Newtonian impact predictions, where Mf = 91.2. It 
is evident from these results that neglecting aerodynamic heating results in large errors and a unrealistic 
prediction of the flutter boundary. With a s = 2°, comparing the heated wing, where Mf = 23.8, to the cold 
wing, where Mf = 62.6, the flutter boundary is reduced by 62% when aerodynamic heating is included. Also, 
the flutter boundary of the heated wing is 75% lower than the flutter boundary of the cold wing predicted 
using Newtonian impact aerodynamics. Based on these limited results it is evident that aerodynamic heating 
is an important consideration, however, this more complicated problem requires additional study. 
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D. Aeroelastic Behavior of a Generic Reusable Launch Vehicle 

The model employed in this study is a vehicle that resembles a generic reusable launch vehicle. It consists 
of a lifting body and canted fins, and has some resemblance to the X-33 vehicle shown in Fig. V. The 
dimensions of the generic vehicle are 76.2 ft. length, 45.54 ft. width, and 6 ft. thickness. The canted fins 
have a span of 18 ft. with a taper ratio of 0.25. They have double-wedge cross-sections with the maximum 
thickness at mid-chord, equal to 3.33% of the chord. The empty mass of the vehicle is considered to be 
70,000 lbs. The unrestrained modes were obtained using MSC.NASTRAN. The first five modes are depicted 
in Fig. 36. The first two modes are the symmetric and antisymmetric fin bending modes, while the third and 
fourth modes show fuselage bending and fuselage torsion. The fifth and higher modes show a combination 
of fin and fuselage deformation. 

In previous hypersonic studies, the approximate aerodynamic method of choice for bluff bodies has been 
Newtonian impact theory blended with piston theory. 16,68 The generic vehicle, shown in Fig. V, has sharp 
leading edges. For M > 5, the wedge angles are small enough to allow the presence of attached oblique 
shocks. However, the wedge angles are large and Mr 1, hence, instead of using a blended approach, a 
different approach, described next, is used. 

The geometry of the generic hypersonic vehicle usually consists of a combination of flat surfaces having 
different wedge angles to the oncoming flow. In regions around the nose, and at the rear of the vehicle, the 
actual slope might be different due to localized curvature of the geometry. Since piston theory is a two- 
dimensional theory which relates parameters in the x-z plane, a sectional model of the vehicle can be created 
using piston theory alone. Fig. 37 shows that the streamlines along the surface in an inviscid flow have 
relatively small flow velocity components in the spanwise direction. Furthermore, at these Mach numbers, 
only a small portion at the root of the fin is affected by the flow around the fuselage (body). Therefore, the 
flow over the fins is assumed to be free-stream in this model. Thus, the surface of the vehicle can be divided 
into “flow zones”, where the flow variables have substantially different values, due to the presence of oblique 
shocks and expansion fans. The local flow variables can thus be calculated at various points on the surface 
of the rigid geometry. Assuming small perturbations during the unsteady aeroelastic analysis, these values 
need not be recalculated, saving computational effort. Next, an unsteady aeroelastic analysis is carried out 
using piston theory with these local variables, applied in frames rotated to the angles of the flat surfaces 
approximating the fuselage. 

Using third order piston theory, the aeroelastic behavior of the generic hypersonic vehicle was studied 
at different altitudes. The flutter Mach number, Mf, was calculated at 40,000 feet using piston theory with 
a varying number of structural modes. The results given in Table 6 indicate that the minimum number of 
modes required for a converged solution is 12. Results obtained at 40,000 ft using piston theory, with 12 
structural modes, are shown in Fig. 38, where Mf =7.1 and fc w =0.3106. In this case, flutter arises from the 
interaction among the higher modes, which have significant fin bending. The fuselage flexural and torsional 
modes (modes 3 and 4) do not play a significant role in this instability. Note, that tracking a substantial 
number of modes required a new graphical approach, 43 as shown in Fig. 38. The upper half of the figure 
shows the variation in damping with Mach number, for the modes which are identified by a mode number 
that is specified along the third axis. The flutter Mach number is predicted when any of the damping curves 
rises above the shaded plane, which corresponds to zero damping. The lower half of the figure shows the 
variation of modal frequencies with Mach number for the various modes. The frequencies of the modes at 
flutter are indicated by black dots. 

The values of the reduced frequencies for this case indicates that the use of a quasisteady approach such 
as piston theory is justified. Note that the reduced frequencies are calculated using the vehicle length as the 
reference length, and would be much smaller if the root chord length of the fin were used. 

Figures 39 and 40 present the aeroelastic behavior of the generic vehicle using Euler aerodynamics. 
Figure 39 shows the aeroelastic behavior of the vehicle calculated on the coarse mesh, while Fig. 40 shows 
the aeroelastic behavior calculated on the refined mesh. Both series of simulations predict flutter at Mf= 9.3, 
which is 31% higher than that predicted by nonlinear piston theory. The figures also show that all modes 
become equally unstable as the flutter boundary is reached, though the higher modes with significant fin 
bending are slightly less stable. The reduced frequency at flutter obtained on the coarse mesh was £^=0.1714, 
while that on the refined mesh was fc w =0.1626. 

The significant difference in flutter boundaries predicted using the two aerodynamic models is due to the 
presence of three-dimensional flow effects. Prior to calculation of the aeroelastic transients, the aeroelastic 
equilibrium state of the flexible vehicle was calculated by allowing the vehicle to deform in hypersonic flow, 


17 of 47 


American Institute of Aeronautics and Astronautics 


with artificial structural damping added to accelerate convergence. In the case of piston theory which is a 
two-dimensional theory, as each cross-section of the vehicle is symmetrical about the horizontal plane, the 
equilibrium state of the flexible vehicle is the same as the rigid geometry at zero angle-of-attack. However, 
Euler aerodynamics is a fully three-dimensional aerodynamic model, and the cantedness of the fins leads to 
higher pressure on the upper surface compared to the lower surface, causing the fins to bend outwards along 
the span-wise direction. For inviscid flow with M= 8.0 at 40,000 feet, the typical cross-section translates 
downward by nearly 1 foot. This deformation of the section is strongly dependent on the specific structural 
model of the fins and stiffeners. 

Figure 41 depicts the typical cross-sections of a fin from rigid and aeroelastic equilibrium states in Euler 
flow, now superimposed. It is seen that the cross-section of the flexible vehicle has deformed and also 
developed a slight angle of attack. Carrying out a linear fit to the mean line of the cross-section in Euler 
flow indicates an effective angle of attack of 0.3°. However, Ref. 25 has indicated that introducing a small 
constant angle of attack in hypersonic flow, for the case of a two-dimensional double-wedge airfoil, produces 
only a slight change in the flutter boundary of the airfoil section. Figure 42 shows the surface pressures 
on the typical cross-sections of the rigid and flexible configurations at equilibrium in inviscid flow, and it is 
evident from the figure that these are quite different. Note that the surface pressures on the rigid vehicle 
using nonlinear piston theory is the same as that obtained using Euler aerodynamics. The deformation of 
the typical cross-section from the double-wedge airfoil leads to a significantly different pressure profile of 
the typical section. Thus, the introduction of flexibility interacts with the three-dimensional flow captured 
using Euler aerodynamics to change the equilibrium configuration of the generic vehicle. In comparison, the 
aeroelastic equilibrium state for the vehicle is the same as the undeformed state when using piston theory. 
Since the aeroelastic transients are calculated by considering small perturbations about the equilibrium 
position, the two models yield different aeroelastic behavior. 

Figure 43 depicts the contours of C v around the typical section for the rigid vehicle at M= 8.6 and 40,000 
feet when using Euler aerodynamics, and Fig. 44 shows the flow around the deformed section when vehicle 
flexibility is introduced. Comparing the two figures shows the significant change in C p on and around the 
typical section. Similar differences were seen at M= 2.6 at 10,000 feet and M— 25.0 at 70,000 feet. As can be 
seen from these figures, the deformed section leads to an asymmetric distribution of pressure on the top and 
bottom surfaces, especially near the leading edge, unlike the rigid section which develops symmetric pressure 
loading on both surfaces. 

For comparison, the aeroelastic simulations of the entire vehicle were also carried out using 7 modes. For 
this incomplete representation of the generic vehicle, both piston-theory and Euler-based aeroelastic studies 
indicate similar flutter behavior, with flutter arising from interaction between the symmetric fin bending 
mode and the fuselage flexural mode. This trend of good comparison between Euler aerodynamics and 
piston theory follows that also observed for the low aspect-ratio wing, and Ref. 19 for the X-43 geometry. 
For convenient comparison, the results of all the aeroelastic simulations described above are summarized in 
Table 7. 

The flutter envelope for the generic vehicle, calculated using 12 modes, is shown in Fig. 45. This 
figure compares the flutter boundaries at various altitudes calculated using third-order piston theory, and 
Euler aerodynamics on the coarse and refined meshes. At low altitudes, nonlinear piston theory predicts 
flutter boundaries comparable with Euler aerodynamics-based methods for the generic vehicle. However, at 
altitudes above 30,000 feet, there is a significant difference, which becomes approximately constant above 
50,000 feet. For the range between 10,000 and 50,000 feet, for which results could be obtained, the coarse 
mesh also compares well with the flutter results from the refined mesh. 

The differences displayed in Fig. 45 are magnified in Fig. 46, which compares the dynamic pressure at 
flutter for the different altitudes considered. From this figure, it is obvious that the comparison between piston 
theory and Euler aerodynamics at low altitudes is not as good as seen in Fig. 45. At higher altitudes, piston 
theory underpredicts the flutter dynamic pressure compared to Euler aerodynamics by 25-50%. However, 
piston theory does indicate a similar trend, including the maximum value corresponding to 70,000 feet. 
Considering the Euler results on the two meshes, good comparison is observed between 10,000 and 40,000 
feet, while the results at 50,000 feet begin to show a divergence. Results at higher altitudes could not be 
obtained on the coarse mesh due to the insufficiently refined grid spacing near the surface, leading to a lack 
of convergence in the flow solutions. Solutions above 80,000 feet could not be obtained on the refined mesh 
due to a similar lack of convergence for M > 28. 

Another challenge with calculating flutter boundaries at very high altitudes is the length of the transient 
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required to identify flutter. At altitudes below 70,000 feet, transients of 1 second have been used to identify 
flutter reliably with the MBA. At very high altitudes, due to the rarefied atmosphere, the aeroelastic loads 
are quite small even at high M. Figure 47 shows the aeroelastic transients for modes 1 and 9, calculated 
at M= 28.0 and 80,000 feet. The very gradual increase in the amplitude of mode 9 can be seen, although 
all the other modes were observed to be quite stable. For this flight condition, the computational resources 
required are an order of magnitude higher than that required for lower altitudes, as the transients have to be 
simulated over much longer periods. This is an excellent example of the importance of choosing an efficient 
time domain damping and frequency identification technique, such as the ARMA method. 

The simulations in this study were carried out in parallel mode on nodes of a AMD cluster with 2000MP 
processors. For the time-step size selected, the run time to obtain an aeroelastic transient response of the 
generic vehicle at 40,000 feet using Euler aerodynamics was approximately 30 hours with 5 processors on the 
coarse mesh, and approximately 140 hours with 9 processors on the refined mesh. The use of either 7 or 12 
structural modes did not significantly affect run time. A sample calculation carried out using Navier-Stokes 
aerodynamics required 36 days on 9 processors when using the refined mesh. 

IV. Conclusions 

The fundamental studies in hypersonic aeroelasticity and aerothermoelasticity allow one to reach several 
useful conclusions. 

1. In general, the three time domain frequency and damping identification techniques produce similar 
estimates of the aeroelastic behavior of a system. The ARMA method was superior, however, to 
both the LSCFM and MBA method since it provided quick damping and frequency estimates with 
minimal response record length. This is an important consideration because the cost of computational 
aeroelastic simulations is largely determined by the length of response record required to accurately 
calculate system damping. In this case, the ARMA method offers a 75% reduction in computational 
cost over the MBA. 

2. The aeroelastic behavior of a system, predicted using time-accurate CFD solutions of the Euler equa- 
tions are sensitive to the number of sub-iterations used, as well as CFL T . Studies must be carried out 
on a case by case basis in order to ensure that these parameters have appropriate values. 

3. The aeroelastic behavior of the three-dimensional low aspect ratio wing obtained using piston theory, 
Euler and Navier-Stokes aerodynamics is similar. Euler solutions are approximately 2% higher than 
those predicted by piston theory, while Navier-Stokes solutions are less than 1% lower than the Euler 
solutions. 

4. The presence of aerodynamic heating on a low aspect ratio cantilever structure, such as a fin and/or 
control surface on a hypersonic vehicle, will result in thermal stresses due to warping restraint at the 
root. This, combined with material property degradation, dramatically affects the natural frequencies 
of the structure. Increasing the Mach number decreases the time before thermal buckling occurs. 
However, changes in frequency are similar for various Mach numbers when plotted as a function of 
leading edge temperature. 

5. Angle of attack is an important consideration when performing aerothermoelastic analysis of a struc- 
ture, since it introduces additional thermal stresses which significantly degrade the stiffness of the 
structure for a given reference temperature, or point in time. Aerodynamic heating in this case has a 
substantial effect on the aeroelastic behavior of the low aspect ratio wing, reducing the flutter boundary 
by 62%. 

6. For the generic vehicle, the aeroelastic model based on Euler solutions predicts a flutter boundary 
31% higher than that predicted by third-order piston theory, due to the presence of significant three- 
dimensional flow effects captured by Euler aerodynamics. This leads to significant deformation of 
the canted fins, and different equilibrium configurations for the flexible vehicle are predicted by these 
two aerodynamic models. Hence, the surface pressure distributions are also found to be different. 
Since the aeroelastic transient solutions are calculated by considering small perturbations about these 
equilibrium positions, different aeroelastic behavior is predicted by these models. 
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7. Based on the various aeroelastic simulations described in this study, it can be seen that the differences 
between approximate solutions based on third-order piston theory and more exact solutions based on 
Euler and Navier-Stokes aerodynamics can be small, or large, and these differences depend on the 
vehicle configuration and flight conditions. 

8. The relatively low reduced frequencies observed for both the low aspect ratio wing and the generic 
hypersonic vehicle can be used to justify the quasisteady nature of piston theory loads. Also, such 
an approximation can also be used to increase computational efficiency when carrying out Euler or 
Navier-Stokes computations. 

9. The aeroelastic behavior of the configurations considered is relatively insensitive to grid resolution as 
long as a sufficient number of grid points are within the shock layer of the flow. 

10. The results presented can be considered to provide a partial validation of the aeroelastic capabilities 
of the CFL3D code for the hypersonic flow regime. 
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Table 1. Error in moment calculation for typical section as a function of grid spacing. 


Nondimensional grid-spacing 
in normal direction (A*.) 

C M 

Error (%) 

2 x 10" 6 

-0.1281 

0.2799 

5 x 10" 7 

-0.1279 

0.1347 

1 x 1CT 8 

-0.1278 

0.0490 

1 x 10" 9 

-0.1277 

N/A 


Table 2. Parameters describing the double- wedge airfoil. 


Parameter 


c (m) 

2.35 

Thickness ratio (%) 

3.36 

Wedge angle (°) 

3.85 

m (kg/m) 

94.2 

r a 

0.484 

u>h (Hz) 

13.4 

u> a (Hz) 

37.6 

Uh 

0.356 

^OL 


•&CX. 

0.2 

a 

0.1 


Table 3. Comparison of the Lockheed F-104 Starfighter wing to the low aspect ratio wing model 


Parameter 

F-104 

3-D Wing 

Wing Mass (Kg) 

350.28 

350.05 

1st Bending Frequency (Hz) 

13.40 

13.41 

1st Torsional Frequency (Hz) 

37.60 

37.51 


Table 4. Flutter results for the low aspect ratio wing using an increasing number of modes, and third order 
piston theory aerodynamics, at 40,000ft. 


No. of Modes 

M f 

2 

15.6 

5 

13.4 

8 

13.3 
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Table 5. Comparison of the original and modified low aspect ratio wings 


Parameter 

Original Wing 

Modified Wing 

Wing Mass (Kg) 

350.05 

377.73 

Mode 1 Freq. (Hz) 

13.41 

14.28 

Mode 2 Freq. (Hz) 

37.51 

40.94 

Mode 3 Freq. (Hz) 

49.18 

60.06 

Mode 4 Freq. (Hz) 

77.14 

81.86 

Mode 5 Freq. (Hz) 

79.48 

97.25 


Table 6. Flutter results for the generic vehicle using an increasing number of modes, and third order piston 
theory aerodynamics, at 40,000ft. 


No. of modes 

M f 

7 

11.6 

10 

7.8 

12 

7.1 

14 

7.1 


Table 7. Comparison of flutter boundaries for the generic vehicle at 40,000 feet, obtained using different 
models for predicting aerodynamic loads. 


Aerodynamic 

model 

7 modes 

12 modes 

M f 

kuj 

M f 

kuj 

Piston theory, nonlinear 

11.7 

0.1344 

7.1 

0.3106 

Euler, coarse mesh 

11.6 

0.1412 

9.3 

0.1714 

Euler, refined mesh 

11.2 

0.1426 

9.3 

0.1626 
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300,000 



Figure 1. Operating envelopes for several modern hypersonic vehicles. 



transient 

response 


Postprocessing 
( , , flutter boundary) 


Figure 2. A flow diagram of the computational aeroelastic solution procedure. 


25 of 47 


American Institute of Aeronautics and Astronautics 



y (m) 


Figure 3. Computational domain of the low aspect ratio wing. 



Figure 4. Span section plane of the computational domain at the root of the low aspect ratio wing. Note that 
the z dimension is scaled relative to the x and y dimensions for visualization purposes. 
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(a) X33. top view. 


(b) Generic vehicle, top view. 



(c) X33. front view. 



(d) Generic vehicle, front view. 




(e) X33, side view. 


(f) Generic vehicle, side view. 


Figure 5. X-33 and generic reusable launch vehicle. 



Figure 6. CFD domain for aeroelastic analysis of the generic hypersonic vehicle. 
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Figure 7. A view of the midplane of the grid around the generic vehicle. 



Figure 8. A view of the grid around the typical cross-section of the generic vehicle. 



Surface / TPS 


Figure 9. Heat transfer at the wall of a hypersonic vehicle. 
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transient 

response 


Postprocessing 
(^, oo, flutter boundary) 


Figure 10. A flow diagram for the computational aerothermoelastic solution procedure. The aerodynamic 
heating information is passed from the CFD solver to the finite element analysis to compute the heat transfer 
between the fluid and structure. 



modes, 

frequency 

\ 


Figure 11. A flow diagram of the finite element analysis procedure when the heat transfer between the fluid 
and structure is included in the analysis. 
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z 



Figure 12. Two degree-of- freedom typical airfoil geometry. 



Mode 1 - ARMA 
Mode 2 - ARMA 
Mode 1 - MBA 
Mode 2 - MBA 
Mode 1 - MBA2 



10 11 

Mach Number 


Figure 13. Aeroelastic behavior, as represented by damping and frequency, for a typical section using Navier- 
Stokes aerodynamics at 40,000 ft; predicted using different time domain identification techniques. 
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x 10" 5 



x 10" 5 



Figure 14. Frequency spectrum of the aeroelastic transients of a typical section at M = 10.0 and 40,000ft. 




Figure 15. Frequency spectrum of the aeroelastic transients of a typical section at M = 11.0 and 40,000ft. 
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1 



Figure 16. Frequency spectrum of the pitch response of a typical section at M = 11.0 and 40,000ft. 



Figure IT. Frequency spectrum of the pitch response of a typical section at M = 11.75 and 40,000ft. 
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Figure 18. Aeroelastic behavior, of a typical section using Navier-Stokes aerodynamics, predicted using an 
ARM A model of the aeroelastic system and different time record lengths of the response, at 40,000 ft. 



Figure 19. Aeroelastic behavior, of a typical section using Navier-Stokes aerodynamics, predicted using MBA 
and different time record lengths of the response, at 40,000 ft. 
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Figure 20. Aeroelastic behavior, of a typical section using Navier-Stokes aerodynamics, predicted using MBA2 
and different time record lengths of the response, at 40,000 ft. 



Figure 21. Aeroelastic behavior, of a typical section using Navier-Stokes aerodynamics, predicted using LSCFM 
and different time record lengths of the response, at 40,000 ft. 
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Flow Direction 



Figure 22. Planform and cross-sectional views of the low aspect ratio wing. 


Mode 1, 13.41 Hz.. First Bending. 



Mode2, 37.51 Hz.. First Torsional. 



Mode 3, 49.18 Hz.. Second Bending. 



Mode 4, 77. 14 Hz.. Second Torsional. 


Mode 5, 79.48 Hz.. Third Torsional. 




Figure 23. First 5 free vibration modes of the low aspect ratio wing 
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Figure 24. Aeroelastic behavior of the low aspect ratio wing using 3rd order piston theory aerodynamics at 
40,000 ft. 



Mode 1 
Mode 2 
Mode 3 
Mode 4 
Mode 5 


CFL t = 1.0 
CFL t = 5.0 


Figure 25. Effect of increasing sub-iterations and CFL T on the aeroelastic behavior of the low aspect ratio 
wing using Euler aerodynamics at 40,000 ft. Symbols indicate the mode number, while the solid vs. dashed 
lines indicate the value of CFL T . 
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Figure 26 . Aeroelastic behavior of the low aspect ratio wing using Euler aerodynamics at 40,000 ft. The 
results were generated using a 0.63xl0 6 cell grid. 



Figure 27 . Aeroelastic behavior of the low aspect ratio wing using Euler aerodynamics at 40,000 ft. The 
results were generated using a 0.27xl0 6 cell grid. 
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Figure 28. Aeroelastic behavior of the low aspect ratio wing using Navier-Stokes aerodynamics at 40,000 ft 

CFLr = 1.0 



Figure 29. Aeroelastic behavior of the low aspect ratio wing using Navier-Stokes aerodynamics at 40,000 ft 
CFLr = 5.0 
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Mode 1, 14.28 Hz. First Bending. Mode2, 40.95 Hz. First Torsional. 



Figure 30. First 5 free vibration modes of the modified low aspect ratio wing 




M = 5.0, Re = 1.3e8 
(T R ) maI = 1065.6 K 


a) 40,000 ft 


M = 10.0, Re = 2.6e8 
(T„)— = 2116.6 K 




M = 5.0, Re = 1.8e7 
(TJma, = 927.4 K 


M = 10.0, Re = 3.7e7 
(TO— = 1650.1 K 




M = 5.0, Re = 6.5e6 M = 10.0, Re = 1.3e7 

(T R ) max = 938.0 K c) 100,000 ft 00- = 1492.7 K 


Figure 31. Normalized radiation equilibrium wall temperatures of the modified low aspect ratio wing for 
various altitudes and free stream Mach numbers. 
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Figure 32. Natural frequencies of the modified low-aspect ratio wing subject to aerodynamic heating at 100,000 
ft. 



Figure 33. Natural frequencies of the modified low-aspect ratio wing subject to aerodynamic heating at 100,000 
ft. Frequencies are plotted as a function of T Ref, which corresponds to the leading edge temperature of the 
75% span location. 



Figure 34. Natural frequencies of the modified low-aspect ratio wing subject to aerodynamic heating, and 
varying angle of attack at 100,000 ft. Frequencies are plotted as a function of Tft e f, which corresponds to the 
leading edge temperature of the 75% span location. 
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Mach Number 


Figure 35. Comparison of flutter boundaries for the modified low aspect ratio wing for the ascent trajectories 
of typical hypersonic aircraft. The flutter boundary for the heated wing at a s =2° corresponds to T Ref = 665 
K in Figure 34. 



First Fuselage Bending 



Mode 4: 14.59 Hz. Mode 5: 21.95 Hz. 

First Fuselage Torsion Second Fuselage Bending 




Figure 36. First 5 free vibration modes of the generic reusable launch vehicle.. 
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Figure 37. Streamlines and Mach contours on the generic vehicle using Euler aerodynamics, top view, M=10.5, 
40,000 feet. 




Figure 38. Aeroelastic behavior of the generic vehicle at 40,000 feet, calculated using nonlinear piston theory, 
and 12 modes. 
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Figure 39. Aeroelastic behavior of the generic vehicle at 40,000 feet, calculated using Euler aerodynamics, 
and 12 modes. The calculations were carried out on the coarse mesh. 






Figure 40. Aeroelastic behavior of the generic vehicle at 40,000 feet, calculated using Euler aerodynamics, 
and 12 modes. The calculations were carried out on the fine mesh. 
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Typical section, rigid structure 


Typical section, flexible structure, superimposed 


□ Mean line, rigid structure 


o Mean line, flexible structure, superimposed 



Figure 41. Superimposed views of the typical sections and mean lines for the rigid and damped equilibrium 
configurations of the typical cross-section. M = 8.0, 40,000 feet. 


a Top surface, rigid structure 

v Bottom surface, rigid structure 

a Top surface, flexible structure, superimposed 

v- Bottom surface, flexible structure, superimposed 



x-coordinate (feet) 


Figure 42. C v on the surface for the rigid and damped equilibrium configurations of the typical cross-section. 
M = 8.0, 40,000 feet. 
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Figure 43. Pressure contours around the rigid typical section of the canted fin of the generic vehicle at M=8.0, 
at 40,000 feet, using Euler aerodynamics. 



Figure 44. Pressure contours around the deformed typical section of the canted fin of the generic vehicle at 
M=8.0, at 40,000 feet, using Euler aerodynamics. 
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Figure 45. Flutter envelope of the generic hypersonic vehicle, calculated using piston theory and Euler aero- 
dynamics. 


— e— 

— 3 rd order piston theory 



• Euler, coarse 

— v- 
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Figure 46. Dynamic pressure at flutter for the generic vehicle, calculated using piston theory and Euler 
aerodynamics. 
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Figure 47. Aeroelastic transients for the generic vehicle at M= 28.0 and 80,000 feet. Only modes 1 and 9 are 
shown for clarity. 
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